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Abstract 

A  new  collocated  finite  volume-based  solution  procedure  for  predicting  viscous 
compressible  and  incompressible  flows  is  presented.  The  technique  is  equally 
applicable  in  the  subsonic,  transonic,  and  supersonic  regimes.  Pressure  is 
selected  as  a  dependent  variable  in  preference  to  density  because  changes  in 
pressure  are  significant  at  all  speeds  as  opposed  to  variations  in  density  which 
become  very  small  at  low  Mach  numbers.  The  newly  developed  algorithm  has 
two  new  features:  (i)  the  use  of  the  Normalized  Variable  and  Space  Formulation 
methodology  to  bound  the  convective  fluxes;  and  (ii)  the  use  of  a  high- 
resolution  scheme  in  calculating  interface  density  values  to  enhance  the  shock 
capturing  property  of  the  algorithm.  The  virtues  of  the  newly  developed  method 
are  demonstrated  by  solving  a  wide  range  of  flows  spanning  the  subsonic, 
transonic,  and  supersonic  spectrum.  Results  obtained  indicate  higher  accuracy 
when  calculating  interface  density  values  using  a  High-Resolution  scheme. 


Nomenclature 


a a Coefficients  in  the  discretized  equation. 

b*  Source  term  in  the  discretized  equation  for  f 

Cp  Coefficient  equals  to  1/RT. 

D[(j)]  The  D  operator. 

df  Covariant  unit  vector  (i.e.  in  the  direction  of  df ). 

D[(t)]  The  vector  form  of  the  D  operator. 

Ff  Convective  flux  at  cell  face  'f . 
ff  Interpolation  factor. 

H[(t)]  The  H  operator. 

H[(j)]  The  vector  form  of  the  H  operator, 
i  Unit  vector  in  the  x-direction. 

j  Unit  vector  in  the  y-direction. 

Jf  Total  scalar  flux  across  cell  face  'f  due  to  convection. 

J°  Total  scalar  flux  across  cell  face  'f  due  to  diffusion. 

Jj.  Total  scalar  flux  across  cell  face  'f . 

M  Mach  number 

P  Pressure. 

Q*  Source  term  for  (j). 

R  Gas  constant. 

Sf  Surface  vector. 

Sf  Contravariant  unit  vector  (i.e.  in  the  direction  of  Sf ). 
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T  Temperature, 

t  Time. 

u,  V  Velocity  components  in  the  x-  and  y-  directions. 
Uf  Interface  flux  velocity  (vf.Sf). 

V  ul  +  vj . 

X,  y  Cartesian  coordinates. 

||a,b||  The  maximum  of  a  and  b. 

Greek  Symbols 


gf 

A[  ([)] 
O 

p<t> 

Q 

a 

P 

5t 

<1) 

Y 

P 

V 


Space  vector  equal  to  (uf  -yd^  )Sf 
The  A  operator. 

Dissipation  term  in  energy  equation. 

Diffusion  for  (j). 

Cell  volume. 

Under-relaxation  factor. 

Thermal  expansion  coefficient  also  flow  angle  at  inlet. 
Time  step. 

Normalized  scalar  variable. 

Scalar  variable. 

Scaling  factor  also  specific  heat  ratio. 

Viscosity. 

Density. 


dx 


d  .  d  I 
—  J  =  — g 
dy 
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Subscripts 


e,  w,  . 

Refers  to 

E,W,.. 

Refers  to 

f 

Refers  to 

F 

Refers  to 

NB 

Refers  to 

P 

Refers  to 

Superscripts 


0 

Refers  to 

(n) 

Refers  to 

* 

Refers  to 

i 

Refers  to 

C 

Refers  to 

D 

Refers  to 

DC 

Refers  to 

DN 

Refers  to 

HR 

Refers  to 

U 

Refers  to 

X 

Refers  to 

y 

Refers  to 

Refers  to 

the  east,  west, ...  face  of  a  control  volume. 

the  East,  West, ...  neighbors  of  the  main  grid  point. 

control  volume  face  f. 

main  grid  point  F. 

neighbours  of  the  P  grid  point. 

the  P  grid  point. 


values  from  the  previous  time  step, 
value  from  the  previous  iteration, 
values  from  the  previous  iteration, 
correction  field, 
convection  contribution, 
diffusion  contribution, 
cross  diffusion, 
normal  diffusion, 
values  based  on  a  HR  scheme, 
values  based  on  the  UPWIND  scheme, 
component  in  x-direction. 
component  in  y-direction. 
dependent  variable. 


Introduction 


In  Computational  Fluid  Dynamics  (CFD)  a  great  research  effort  has  been 
devoted  to  the  development  of  accurate  and  efficient  numerical  algorithms 
suitable  for  solving  flows  in  the  various  Reynolds  and  Mach  number  regimes. 
The  type  of  convection  scheme  to  be  used  in  a  given  application  depends  on 
the  value  of  Reynolds  number.  For  low  Reynolds  number  flows,  the  central 
difference  or  hybrid  scheme  is  adequate  [1].  In  dealing  with  flows  of  high 
Reynolds  number,  numerous  discretization  schemes  for  the  convection  term 
arising  in  the  transport  equations  have  been  employed  [2,3,4,5,6,7,8,9,10,11]. 
On  the  other  hand,  the  Mach  number  value  dictates  the  type  of  algorithm  to  be 
utilized  in  the  solution  procedure.  These  algorithms  can  be  divided  into  two 
groups:  density-based  methods  and  pressure-based  methods,  with  the  former 
used  for  high  Mach  number  flows,  and  the  latter  for  low  Mach  number  flows.  In 
density-based  methods,  continuity  is  employed  as  an  equation  for  density  and 
pressure  is  obtained  from  an  equation  of  state,  while  in  pressure-based 
methods,  continuity  is  utilized  as  a  constraint  on  velocity  and  is  combined  with 
momentum  to  form  a  Poisson  like  equation  for  pressure.  Each  of  these 
methods  is  appropriate  for  a  specific  range  of  Mach  number  values. 

The  ultimate  goal,  however,  is  to  develop  a  unified  algorithm  capable  of  solving 
flow  problems  in  the  various  Reynolds  and  Mach  number  regimes.  To 
understand  the  difficulty  associated  with  the  design  of  such  an  algorithm,  it  is 
important  to  understand  the  role  of  pressure  in  a  compressible  fluid  flow  [12].  In 
the  low  Mach  number  limit  where  density  becomes  constant,  the  role  of 
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pressure  is  to  act  on  velocity  through  continuity  so  that  conservation  of  mass  is 
satisfied.  Obviously,  for  low  speed  flows,  the  pressure  gradient  needed  to  drive 
the  velocities  through  momentum  conservation  is  of  such  magnitude  that  the 
density  is  not  significantly  affected  and  the  flow  can  be  considered  nearly 
incompressible.  Hence,  density  and  pressure  are  very  weakly  related.  As  a 
result,  the  continuity  equation  is  decoupled  from  the  momentum  equations  and 
can  no  longer  be  considered  as  the  equation  for  density.  Rather,  it  acts  as  a 
constraint  on  the  velocity  field.  Thus,  for  a  sequential  solution  of  the  equations, 
it  is  necessary  to  devise  a  mechanism  to  couple  the  continuity  and  momentum 
equations  through  the  pressure  field.  In  the  hypersonic  limit  where  variations  in 
velocity  become  relatively  small  as  compared  to  the  velocity  itself,  the  changes 
in  pressure  do  significantly  affect  density.  In  this  limit,  the  pressure  can  be 
viewed  to  act  on  density  alone  through  the  equation  of  state  so  that  mass 
conservation  is  satisfied  [12]  and  the  continuity  equation  can  be  viewed  as  the 
equation  for  density.  This  view  of  the  two  limiting  cases  of  compressible  flow 
can  be  generalized  in  the  following  manner.  In  compressible  flow  situations,  the 
pressure  takes  on  a  dual  role  to  act  on  both  density  and  velocity  through  the 
equation  of  state  and  momentum  conservation,  respectively,  so  that  mass 
conservation  is  satisfied.  For  a  subsonic  flow,  mass  conservation  is  more 
readily  satisfied  by  pressure  influencing  velocity  than  pressure  influencing 
density.  For  a  supersonic  flow,  mass  conservation  is  more  readily  satisfied  by 
pressure  influencing  density  than  pressure  influencing  velocity. 

The  above  discussion  highlights  the  difficulties  associated  with  the  use  of 
density  as  a  primary  variable  for  computing  low  Mach  number  flows  or  mixed 
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compressible  and  incompressible  flows.  Most  importantly,  it  reveals  that  for  any 
numerical  method  to  be  capable  of  predicting  both  incompressible  and 
compressible  fluid  flows  the  pressure  should  always  be  allowed  to  play  its  dual 
role  and  to  act  on  both  velocity  and  density  to  satisfy  continuity.  Hence, 
pressure  methods  developed  for  incompressible  flow  cannot  be  directly  used 
for  simulating  flow  at  high  Mach  number  without  proper  modifications. 
Furthermore,  the  behavior  of  pressure  in  compressible  flow  depends  on  Mach 
number.  In  the  subsonic  regime,  it  exhibits  an  elliptic  behavior  and  the 
disturbances  at  any  location  affect  and  are  affected  by  all  neighboring  points.  In 
the  supersonic  regime,  however,  the  behavior  is  hyperbolic  and  the  value  at  a 
point  depends  on  the  data  bounded  by  the  characteristics  passing  through  that 
point.  The  discretization  scheme  chosen  for  the  pressure  equation  needs  to 
model  these  features  properly. 

Several  researchers  [12,13,14,15,16,17,18,19,20,21,22,23]  have  worked  on 
extending  the  range  of  pressure-based  methods,  with  various  degrees  of 
success,  to  high  Mach  numbers  following  either  a  staggered  grid  approach  [12- 
15]  or  a  collocated  variable  formulation  [16-23].  The  method  of  Shyy  and  Chen 
[14],  developed  within  a  multigrid  environment,  uses  a  second-order  upwind 
scheme  in  discretizing  the  convective  terms.  Moreover,  at  high  Mach  number 
values,  a  first  order  upwind  scheme  is  employed  for  evaluating  the  density  at 
the  control  volume  faces.  Yang  et  al  [16]  used  a  general  strong  conservation 
formulation  of  the  momentum  equations  that  allows  several  forms  of  the  velocity 
components  to  be  chosen  as  dependent  variables.  In  the  method  developed  by 
March!  and  Maliska  [17],  values  for  density,  convection  fluxes,  and  convection- 
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like  terms  at  the  control  volume  faces  are  calculated  using  the  upwind  scheme. 
Demirdzic  et  al  [18],  however,  used  a  central  difference  scheme  blended  with 
the  upwind  scheme  to  evaluate  these  quantities.  Lien  and  Leschziner  [19,20] 
adopted  the  streamwise-directed  density-retardation  concept,  which  is 
controlled  by  Mach-number-dependent  monitor  functions,  to  account  for  the 
hyperbolic  character  of  the  conservation  laws  in  the  transonic  and  supersonic 
regimes.  Politos  and  Giannakoglou  [21]  developed  a  pressure-based  algorithm 
for  high-speed  turbomachinery  flows  following  also  the  retarded  density 
concept.  In  their  method,  unlike  the  work  of  Lien  and  Leschziner  [19,20],  the 
retarded  density  operates  only  on  the  velocity  component  correction  during  the 
pressure  correction  phase.  Chen  and  Pletcher  [22]  developed  a  coupled 
modified  strongly  implicit  procedure  that  uses  the  strong  conservation  forms  of 
Navier-Stokes  equations  with  primitive  variables.  The  method  of  Karimian  and 
Schneider  [23]  is  formulated  within  a  control-volume  finite  element  framework. 
From  the  aforementioned  literature  review,  it  is  obvious  that  in  most  of  the 
published  work  the  first  order  upwind  scheme  is  used  to  interpolate  for  density 
when  in  the  source  of  the  pressure  correction  equation,  exception  being  in  the 
work  presented  in  [18-21]  where  a  central  difference  method  is  adopted.  In  the 
technique  developed  by  Demirdzic  et  al  [18],  the  second  order  central 
difference  scheme  blended  with  the  upwind  scheme  is  used.  The  bleeding 
relies  on  a  factor  varying  between  0  and  1 ,  which  is  problem  dependent  and 
has  to  be  adjusted  to  eliminate  oscillation  or  to  promote  convergence.  In  the 
work  presented  in  [19-21],  the  retarted  density  concept  is  utilized  in  calculating 
the  density  at  the  control  volume  faces.  This  concept  is  based  on  factors  that 
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are  also  problem  dependent  and  requires  the  addition  of  some  artificial 
dissipation  to  stabilize  the  algorithm  (second-order  terms  were  introduced), 
which  complicate  its  use. 

To  this  end,  the  objective  of  this  paper  is  to  present  a  newly  developed 
pressure-based  solution  procedure  that  is  equally  valid  at  all  Reynolds  and 
Mach  number  values.  The  collocated  variable  algorithm  is  formulated  on  a  non- 
orthogonal  coordinate  system  using  Cartesian  velocity  components.  The 
method  is  easy  to  implement,  highly  accurate,  and  does  not  require  any  explicit 
addition  of  damping  terms  to  stabilize  it  or  to  properly  resolve  shock  waves. 
Moreover,  the  algorithm  will  have  two  new  features.  The  first  one  is  the  use  of 
the  Normalized  Variable  Formulation  (NVF)  [24]  and/or  the  Normalized  Variable 
and  Space  Formulation  (NVSF)  [25]  methodology  in  the  discretization  of  the 
convective  terms.  To  the  authors’  knowledge,  the  NVF/NVSF  methodologies 
have  never  been  used  to  bound  the  convective  flux  in  compressible  flows. 
Mainly  low  order  schemes  or  the  TVD  [26]  formulation  has  usually  been 
adopted.  The  second  one  is  the  use  of  High-Resolution  (HR)  schemes  in  the 
interpolation  of  density  appearing  in  the  mass  fluxes  in  order  to  enhance  the 
shock  capturing  capability  of  the  method. 

In  what  follows  the  governing  equations  for  compressible  flows  are  presented 
and  their  discretization  detailed  so  as  to  lay  the  ground  for  the  derivation  of  the 
pressure-correction  equation.  Then,  a  brief  description  of  the  various  types  of 
boundary  conditions  for  both  incompressible  and  compressible  flows  is  given. 
Finally,  the  increase  in  accuracy  with  the  use  of  HR  schemes  for  density  is 
demonstrated.  This  is  done  by  comparing  predictions,  for  a  number  of 
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problems,  obtained  using  the  third-order  SMART  scheme  [8]  for  all  variables 
except  density  (for  which  the  Upwind  [1]  scheme  is  used)  against  another  set  of 
results  obtained  using  the  SMART  scheme  for  all  variables  including  density. 
Results  generated  using  a  number  of  other  HR  schemes  are  also  presented. 

Governing  Equations 


The  equations  governing  the  flow  of  a  two-dimensional  compressible  fluid  are 
the  continuity  equation,  the  momentum  equations,  and  the  energy  equation. 
This  set  of  non-linear,  coupled  equations  is  solved  for  the  unknowns  p,  v,  T  and 
P.  In  vector  form,  these  equations  may  be  written  as: 


^  +  V.(pv)  =  0 
at 

+  V  •  (pw)  =  -VP  +  V  •  (pVv)  +  -  V(pV  •  v) 


g(pT) 

dt 


+  V.(pvT)  =  —  ]  V  •  (kVT)  +  PT 

S  I 


dt 


+  V-(Pv)-PV-(v) 


+  <D  +  q 


where 


3)  =  iJ, 


dx 


+ 


dv 

dy) 


^du  dv^^ 
ydy  dx 


-f(V.vf 


and  p  the  thermal  expansion  coefficient  which  is  equal  to  1/T  for  an  ideal  fluid. 
In  addition  to  the  above  differential  equations,  an  auxiliary  equation  of  state 
relating  density  to  pressure  and  temperature  (p=f(P,T))  is  needed.  For  an  ideal 
fluid,  this  equation  is  given  by: 
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where  R  is  the  gas  constant. 

A  review  of  the  above  differential  equations  reveals  that  they  are  similar  in 
structure.  If  a  typical  representative  variable  is  denoted  by  (j),  the  general 
differential  equation  may  be  written  as, 

^  +  V.(py<,)  =  V.(r*V(|,)+Q‘ 

where  the  expressions  for  r'*’and  Q'*’  can  be  deduced  from  the  parent 
equations.  The  four  terms  in  the  above  equation  describe  successively 
unsteadiness,  convection  (or  advection),  diffusion,  and  generation/dissipation 
effects.  In  fact,  all  terms  not  explicitly  accounted  for  in  the  first  three  terms  are 
included  in  the  catchall  source  term  O’**. 

Finite  Volume  Discretization 

The  general  transport  equation  (Eq.  (6))  is  discretized  using  the  control  volume 
methodology.  The  discretization  process  is  a  two-level  procedure.  In  level  I,  the 
equations  are  integrated  over  a  control  volume  so  as  to  get  a  discretized 
description  of  the  conservation  laws.  In  level  II,  an  interpolation  profile  is  used 
to  relate  the  discretized  terms  of  level  I  to  the  discrete  nodes  in  the  solution 
domain. 

Level  I  discretization 

The  integral  form  of  equation  (6),  obtained  by  applying  the  divergence  theorem 
over  the  control  volume  P  shown  in  Fig.  1 ,  may  be  written  as: 
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^(pv(|)-r*V(t))-ds=  JqMq 

s  n 


where  Q  is  the  volume  of  cell  P.  Replacing  the  surface  integral  over  the  control 
volume  by  a  discrete  summation  of  the  flux  terms  over  the  sides  of  the  control 
volume,  equation  (7)  becomes: 


^[(p4]n+A[(pv.t-r*v<,).s],  =Q‘n 

Ol 

In  the  above  equation,  the  discretized  form  of  the  unsteady  term  will  be  detailed 
later,  and  the  operator  is  the  discretized  version  of  the  surface  integral 
defined  by: 

=<l’e  +<I>w  +<I>s 

Hence  equation  (8)  can  be  written  as 

|.[(p(|.)p]£5  + a[j],  =|-[(p<ti),lQ  +  (J.  +J,  +J.  +  = 

at  at 

where  J  f  represents  the  total  flux  of  (j)  across  face  'f  and  is  given  by 

j,  =  (pv(|)-r‘v4'S, 

The  flux  Jf  is  a  combination  of  the  convection  flux  =  (pv(j))f.Sf  and  diffusion  flux 
jD=(.r<t-V(t))f.Sf. 


Level  11  discretization 

From  equation  (10),  it  is  obvious  that  the  total  fluxes  are  needed  at  the  control 
volume  faces  where  the  values  of  the  dependent  variables  are  not  available 
and  should  be  obtained  by  interpolation.  Therefore,  the  accuracy  of  the  solution 
depends  on  the  proper  estimation  of  these  values  as  a  function  of  the 
neighboring  ^  node  values.  Through  the  use  of  an  interpolation  profile,  or  an 
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estimate  of  how  ^  varies  between  nodes,  the  approximation  scheme  produces 
an  expression  for  the  face  value  which  is  dependent  on  the  nodal  cj)  values  in 
the  vicinity  of  the  face.  Details  regarding  the  profile  assumption  are  presented 
next. 

Discretization  of  the  Unsteady  Term 


For  the  representation  of  the  unsteady  term,  the  grid-point  value  of  (j)  is 
assumed  to  prevail  throughout  the  control  volume.  Then,  as  used  earlier. 


r^dn  =  -®[(p4]n 


The  time  derivative  is  approximated  using  the  following  Euler-implicit 
formulation: 

at  ot 

In  the  above  equation,  6t  represents  the  time  step  and  the  superscript  o 
denotes  values  at  time  (t-5t). 


Discretization  of  the  Diffusion  Fiux  Jf 

The  diffusion  flux  Jf  is  discretized  along  each  surface  of  the  control  volume 
using  the  method  described  in  Zwart  et  al.  [27]  according  to  which  it  is 
decomposed  into: 

(-r‘v+),.s,  =(-r‘V(|.),.n,s,  =-r*[(vrt.(yi),  +(^),.(n,  -(rd),)k 
where  (vrlrjfis  the  average  of  the  adjacent  cell  gradients,  fifand  (Fig.  2)  are 
the  contravariant  (surface  vector)  and  covariant  (curvilinear  coordinate)  unit 
vectors  respectively  ,  and  y  is  a  scaling  factor.  This  factor  is  chosen  such  that  it 
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is  equal  to  1  on  orthogonal  meshes  in  order  for  the  method  to  collapse  to 
classical  stencils  [27,28].  With  that  constraint,  the  expression  for  y  on  structured 
meshes  is: 


1  _  Sf  df 

Il|-  ad^  Sj-  adf 


(15) 


Defining  the  space  vector  e  ^  as: 

s,=(4,-(rd),)s,  =  K;i+K>'j  (16) 

the  expression  for  (-r*V(|i),.S,  becomes, 


(-r*V(|>),.S, 


pij) 

^  f 


(17) 


In  this  form,  the  term  (Vcj))^  .(d)f.  represents  the  gradient  in  the  direction  of  the 
coordinate  line  joining  P  and  F  (see  Figure  2).  Therefore,  the  above  equation 
can  be  rewritten  as: 


(-r‘V(|>),.s,  =-r* 


and  upon  simplifying,  it  reduces  to: 


(18) 


(-r‘vi)),.s,=-r<' 


-<t>p)^^+(v(|))f.(K;ji+Kn) 


S  I*  ^ 


The  (V(t))f  is  calculated  as: 
(^),=f,(V(t))p+(l-f,Xv<t.)p 


(19) 


(20) 


and  the  gradient  at  the  main  grid  point  F  (F=  P,  E,  W,  N,  or  S)  is  obtained  from: 
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(V'l'),  =i|(|«iS  =  i(K.Sp.  +<|,p.S,.  +^>,.Sp.  +4,p.S,.) 

+'I>f,se, 

=(v+)'i+(vrtj 

The  final  form  of  the  diffusive  flux  along  face  f  (T=e,  w,  n,  or  s)  is: 


(-r*vrt.s,  = -r;(+, 


~r*  {jf,  (V'i>);  +  (i  -  f,  Xv'i'):  U  +  if,  (vtf + (1  -  f.  Xv'i’): 


K 


y 

f 


The  underlined  part  of  the  diffusion  flux  is  called  the  cross-diffusion  contribution. 
It  vanishes  when  the  grid  is  orthogonal,  and  is  small  compared  to  normal 
diffusion  for  nearly  orthogonal  grid.  In  such  circumstances,  explicit  treatment  of 
the  cross  diffusion  term  does  not  significantly  influence  the  rate  of  convergence 
of  the  overall  solution  procedure  and  simplifies  the  matrix  of  coefficients. 


Discretization  of  the  Convection  Flux 

The  convection  flux  of  ^  through  the  control  volume  face  f  is  given  by  : 

Jf  =(pv<t))f  -Sf 

where  (t)f  stands  for  the  mean  value  of  (j)  along  cell  face  f,  and  Ff=  (pv.S)f  is  the 
mass  flow  rate  across  face  f.  Using  some  assumed  interpolation  profile,  (t)f  can 
be  explicitly  formulated  by  a  functional  relationship  of  the  form: 

‘t>f  =  f(‘t>„b) 

where  (|)nb  denotes  the  (j)  values  at  the  neighboring  nodes.  The  interpolation 
profile  to  be  used  should  be  bounded  in  order  not  to  give  rise  to  the  well-known 
dispersion  error  problem  [2].  In  this  work,  HR  schemes  formulated  in  the 
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context  of  the  NVSF  methodology,  which  is  explained  in  the  next  section,  are 
used.  Since  these  interpolation  profiles  may  involve  a  large  number  of 
neighboring  points,  the  solution  of  the  resultant  system  of  equations  can 
become  very  expensive  computationally,  hence  the  use  of  a  compacting 
procedure  is  most  welcome.  For  that  purpose,  the  deferred  correction  method 
[29,30]  is  adopted.  With  this  approach,  the  convection  is  split  into  an  implicit 
part,  expressed  through  first  order  upwind  differencing  scheme  (LIDS)  [1],  and 
an  explicit  part,  which  equals  the  difference  between  the  UDS  and  HR 
approximations,  i.e.: 

FA=F,<,“+F,(.|,™-4,“) 

The  deferred  correction  approach  enhances  the  diagonal  dominance  of  the 
matrix  of  coefficient,  which  adds  to  the  stability  of  the  solution  algorithm. 

Discretization  of  the  Sources 

The  integral  value  of  the  source  term  over  the  control  volume  P  (Fig.  1)  is 
obtained  by  assuming  the  estimate  of  the  source  at  the  control  volume  center  to 
represent  the  mean  value  over  the  whole  control  volume.  Hence,  one  can  write: 

JqMQ  =  Q*Q 

n 

Following  the  practice  adopted  in  Patankar  [1],  the  source  term  is  linearized 
according  to: 

Q*  =Q*  +Q*(j) 

where  Q*  should  always  be  a  negative  quantity.  Moreover,  the  additional 
terms,  appearing  in  the  momentum  and  energy  equations,  not  featured  in 
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equation  (1 1 ),  are  treated  explicitly  and  their  discretization  is  analogous  to  that 
of  the  ordinary  diffusion  flux. 

Algebraic  System  Of  Equations 


The  discretized  equation,  Eq.  (10),  is  transformed  into  an  algebraic  equation  at 
the  main  grid  point  P  by  substituting  the  fluxes  at  all  faces  of  the  control  volume 
by  their  equivalent  expressions.  Then,  performing  some  algebraic 
manipulations  on  the  resultant  equation,  the  following  algebraic  relation,  linking 
the  value  of  the  dependent  variable  at  the  control  volume  center  to  the 
neighboring  values,  is  obtained: 

4^,  =  Z^’nb'I'nb  +b;  (28) 

NB(P) 

In  the  above  equation,  are  the  coefficients  multiplying  the  value  of  (j)  at  the 


neighboring  nodes  NB=(E,  W,  N,  and  S)  surrounding  the  central  node  P,  is 
the  coefficient  of  (j)?,  and  b*  contains  all  terms  that  are  not  expressed  through 
the  nodal  values  of  the  dependent  variable  (e.g.  the  source  term  Q'*’,  the 
pressure  gradients  in  the  momentum  equations,  terms  involving  known  values 
of  (|)  etc. ...).  The  final  form  of  these  coefficients  is  as  follows: 


^  -Tf  \~e  /  '  \-e  /  I  II  p 

sM:+srdr  " 


w  w 


-f-  —  p  I 


a’’  =  r''’  _F  0 

“  "sx+sw  ' 


.rtkilkl 


_F  0 

SgXdX^SJdJ  " 


(29) 

(30) 

(31) 


(32) 
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a*  =a* +at +at +a*+a°-Q*Q 

"  At 

b*  =Q*n  +  aJ4,;  -[F.(,|,r  -4,“)+F.(+r  -'t'“)+F.('l'r  -'!>")] 

H-  r‘  {[f.  (V4,);  +  (i  -  f.  Xvf)*  ]<  +  [f.  (vrt + (i  -  f.  Xv<,)’  ]k>  } 

+ rj  {[f.  (v<.X  +  (1  -  f.  XV'I');  ]k:  +  [f,  (V4.)J  +  (1  -  f.  XV'I'fe  Ik^  } 

+r.*  ff,(vrt  +(i-f„Xv+);]K;  +[f.(vrt  +(i-f.Xv^,);]K:} 

+r,‘  ff.(vrt  +(i-f.Xv4,)>]<  +[f,(vrt  +(i-f,Xv4.)']K.>} 

For  the  solution  domain  as  a  whole  there  results  a  system  of  N  equations  in  N 

unknowns,  where  N  is  the  number  of  control  volumes.  Many  techniques  exist 
for  solving  large  systems  of  linear  equations  that  may  be  classified  as  direct  or 
iterative  methods.  The  use  of  direct  methods  is  not  appropriate  in  the  present 
context  because  they  require  much  more  storage  than  iterative  methods  and 
are  usually  more  expensive  computationally.  Owing  to  the  non-linear  nature  of 
the  set  of  equations,  the  discretized  equations  are  solved  by  the  use  of  iterative 
methods.  Current  iterative  methods  differ  with  respect  to  storage  requirement 
and  degree  of  implicitness,  such  as  the  point-by-point  successive  over¬ 
relaxation  method  [31],  the  strongly  implicit  procedure  of  Stone  [32]  and  its 
variations,  the  Incomplete  Cholesky  Conjugent  Gradient  (ICCG)  [33],  or  the 
Multigrid  Method  of  Brandt  [34]  to  site  a  few.  Although  these  methods  have 
their  own  desirable  attributes,  the  degree  of  simplicity  of  their  implementation  in 
a  computer  code  is  approximately  inversely  proportional  to  their  rate  of 
convergence.  The  algorithm  used  in  this  work  is  the  TDMA  [35]. 

During  the  iterative  process,  it  is  often  desirable  to  slow  the  changes,  from 
iteration  to  iteration,  in  the  values  of  the  dependent  variable.  This  process  is 
called  under-relaxation.  It  is  an  important  tool  that  prevents  divergence  of  the 
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iterative  solution  for  strongly  nonlinear  problems,  as  is  the  case  here.  If  (j)*  and 
(j)p  are  the  values  from  the  previous  and  current  iterations,  respectively,  then  Eq. 
(28)  can  be  written  as 


(|)p  —  (|)p  + 


NB(P) 


In  the  above  equation,  the  underlined  term  represents  the  change  in  (j)p 
produced  by  the  current  iteration.  This  change  can  be  varied  by  introducing  an 
under-relaxation  factor  a  (0<a<1 ),  so  that 


(t)p  =  <t)p  +  ot 


^^nb4*nb  ^P 

NB(P) 

i,  Vp 


or 

— <t>P  =  +(l-a)^([); 

a  NBfp^  cc 


At  the  state  of  convergence,  and  (i)p  are  equal  and  the  original  equation  is 
satisfied.  There  is  no  general  rule  for  choosing  the  optimum  value  of  a  and  a 
suitable  one  for  a  given  problem  is  usually  found  from  exploratory 
computations.  Equation  (38)  can  be  rewritten  in  the  form  of  equation  (28)  by 
redefining  a*  and  b*  as  follows: 
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The  NVSF  Methodology  for  Constructing  HR  schemes 

As  mentioned  earlier,  the  discretization  of  the  convection  flux  is  not 
straightforward  and  requires  additional  attention.  Since  the  intention  is  to 
develop  a  high-resolution  algorithm,  the  highly  diffusive  first  order  UPWIND 
scheme  [1]  is  excluded.  As  such,  a  high  order  interpolation  profile  is  sought. 

The  difficulties  associated  with  the  use  of  such  profiles  stem  from  the  conflicting 
requirements  of  accuracy,  stability,  and  boundedness.  Solutions  predicted  with 
high  order  profiles  tend  to  provoke  oscillations  in  the  solution  when  the  local 
Peclet  number  is  high  in  combination  with  steep  gradients  of  the  flow 
properties.  To  suppress  these  oscillations,  many  techniques  have  been 
advertised  and  may  be  broadly  classified  into  two  groups:  the  flux  blending 
method  [36,37,38,39]  and  the  composite  flux  limiter  method  [8,24-26,40],  the 
latter  being  the  one  adopted  here.  In  this  technique,  the  numerical  flux  at  the 
interface  of  the  computational  cell  is  modified  by  employing  a  flux  limiter  that 
enforces  boundedness.  The  formulation  of  high-resolution  flux  limiter  schemes 
on  uniform  grid  has  recently  been  generalized  by  Leonard  [24,40]  through  the 
Normalized  Variable  Formulation  (NVF)  methodology  and  on  non-uniform  grid 
by  Darwish  and  Moukalled  [25]  through  the  Normalized  Variable  and  Space 
Formulation  (NVSF)  methodology.  The  NVF  and  NVSF  methodologies  have 
provided  a  good  framework  for  the  development  of  HR  schemes  that  combine 
simplicity  of  implementation  with  high  accuracy  and  boundedness.  Moreover,  to 
the  authors’  knowledge,  the  NVSF  formulation  has  never  been  used  to  bound 
the  convection  flux  in  compressible  flows.  It  is  an  objective  of  this  work  to 


A  High  Resolution  Pressure  Based  Algorithm  for  Fluid  Flow  at  All  Speeds 


extend  the  applicability  of  this  technique  to  compressible  flows.  Therefore, 
before  introducing  the  high-resolution  algorithm,  a  brief  review  of  the  NVSF 
methodology  is  in  order.  This  is  done  by  first  defining  the  Normalized  Variables, 
then  presenting  the  Convection  Boundedness  Criterion  (CBC)  [8],  and  finally 
describing  the  various  HR  schemes  employed  in  this  work  (OSHER  [41], 
MUSCL  [42],  CLAM  [43],  SMART  [8],  EXPONENTIAL  [24],  SUPER-C  [44],  and 
ISNAS  [45]). 

Normalized  Variables 

Fig.  3(a)  shows  the  local  behavior  of  the  convected  variable  near  a  control- 
volume  face.  The  node  labeling  refers  to  the  upstream,  central,  and 
downstream  grid  points  designated  by  U,,C,  and  D,  located  at  distances  ^u. 
and  4d  from  the  origin,  respectively.  The  values  of  (j)  at  these  nodes  are 
designated  by  (t)u,  (|)c  and  (t)D  respectively.  Moreover,  the  value  of  the  dependent 
variable  at  the  control  volume  face  located  at  a  distance  from  the  origin  is 
expressed  by  (j)f.  With  this  notation,  the  normalized  variables  are  defined  as 
follows: 

The  use  of  the  above-normalized  parameters  simplifies  the  functional 
representation  of  interpolation  schemes  (Fig.  3(b))  and  helps  defining  the 
stability  and  boundedness  conditions  that  they  should  satisfy.  In  addition,  the 
normalized  functional  relationship  of  any  scheme  may  be  plotted  on  a 
Normalized  Variable  Diagram  (NVD)  (Fig.  3(c)),  which  is  an  effective  tool  in 
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assessing  the  accuracy,  boundedness,  and  relative  diffusivity  of  convective 
schemes.  In  general,  the  value  of  (j)f  is  represented  by  the  following  parametric 
relation: 

(|)f  =  f((|)u,(|)c)<|)D»4u’^C>^f  >^d) 
which,  upon  normalizing,  is  simplified  to 

=f($c.5c.5f) 

By  comparing  equations  (41 )  and  (42)  it  is  clear  that  one  of  the  normalization 
benefits  is  to  reduce  the  number  of  parameters  involved  in  the  functional 
relationship. 

The  Convective  Boundedness  Criterion  (CBC) 


Based  on  the  normalized  variable  analysis,  Gaskell  and  Lau  [8]  formulated  a 
convection  boundedness  criterion  (CBC)  for  implicit  steady  flow  calculation. 
This  CBC  states  that  for  a  scheme  to  have  the  boundedness  property  its 
functional  relationship  should  be  continuous,  should  be  bounded  from  below  by 

=  f c  >  above  by  unity,  and  should  pass  through  the  points  (0,0)  and 
(1,1),  in  the  monotonic  range  (0<  fc  for  1<fc  or  <0,  the  functional 


relationship  f(fc  )  should  equal  f  ^  .  Mathematically,  these  conditions  are 


CBC  =  <^ 


m 

m = 0 

f($)=i 

l<f($)<(t)c 


[m = ?c 


is  continuous 
for  ([)(,=  0 
for(j)c  =1 
for  0  <  (j)c  <1 
for  ([)(,  <  0and(t)(,  >  1 


A  High  Resolution  Pressure  Based  Algorithm  for  Fluid  Flow  at  All  Speeds 


24 


The  above  conditions  may  also  be  illustrated  on  a  Normalized  Variable  Diagram 
(NVD)  as  shown  in  Fig.  3(c). 

Normalized  Variable  and  Space  Formulation  (NVSF)  methodology 


Knowing  the  required  conditions  for  boundedness,  the  shortcomings  of  the  High 
Order  (HO)  schemes  were  eliminated  through  the  development  of  HR  schemes 
satisfying  all  above  requirements.  Without  going  into  details,  a  number  of  HR 
schemes  were  formulated  using  the  NVF/NVSF  methodologies  and  the 
functional  relationships  for  the  ones  used  in  this  work  are  given  below.  For 
more  details  the  reader  is  referred  to  Darwish  and  Moukalled  [25]. 


OSHER 


= 


Sc  ^  Sf 

1 


<f'c 


.  -‘t’c  <1 
elsewhere 


(44) 


MUSCL 


<t>c 


?€  +(4f  “^c) 


1 

.‘I’c 


0<?c  <Y 

l  +  ?c-?f 
elsewhere 


(45) 
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CLAM 


5c(?c-l) 


0<(j)c  <1 
elsewhere 


EXPONENTIAL 


j_g-2.19722?cj  0<fc<l 
elsewhere 


SMART 


'ff(l-3fc+2?f) 


<{)f  - 


<t>c 


0  <  <|)c  < 


?c(l-?c) 


$c 


3  •“ 

?c 


“  ^c)  -  ?c 

Sf 


elsewhere 


SUPER-C 


*t>f  - 


<t>c 


^f(l-3|c+25f) 

5c(i“5c) 


i-5c  ■  i-5c 
?f(i-5f).  4f(5f-5c) 


0<(j)c  <  — 


■^C  -  *l>c  < 


<t>C  + - ^ - -  -‘t>c  +  “^c) 

l-4c 


iic(l“Ic) 

1 

l^c 


k 

elsewhere 


“  ^c)  <  ^'c  <  1 

Sf 
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ISNAS 


f  il  -IcJ  r  .  (2?c’  +  -Icl  +?,'  j 


(4c’-34cg,+g/+4,) 


,2  ‘l^C 


0  <  (j)c  <  1 


elsewhere 


High  Resolution  Algorithm 


The  need  for  a  solution  algorithm  arises  in  the  simulation  of  flow  problems 
because  a  scalar  equation  does  not  exist  for  pressure.  Rather,  the  pressure 
field  acts  indirectly  on  the  velocity  field  to  constrain  it  to  satisfy  the  continuity 
equation.  Hence,  if  a  segregated  approach  is  to  be  adopted,  coupling  between 
the  u,  V,  p,  and  P  primitive  variables  in  the  continuity  and  momentum  equations 
will  be  required.  Evidently,  the  whole  set  of  equations  could  be  solved  directly 
(after  linearization)  since  the  number  of  equations  equals  the  number  of 
unknowns.  However,  the  computational  effort  and  storage  requirements  needed 
by  such  an  approach  are  often  prohibitive.  This  has  forced  researchers  to  seek 
less  expensive  methods  and  resulted  in  the  development  of  several  segregated 
solution  algorithms  [1 ,46,47,48,49,50,51 ,52].  Recently,  Moukalled  and  Darwish 
presented  a  unified  formulation  of  these  algorithms  [53]. 

The  segregated  algorithm  adopted  in  this  work  is  the  SIMPLE  algorithm  [1,46], 
that  involves  a  predictor  and  a  corrector  step.  In  the  predictor  step,  the 
velocity  field  is  calculated  based  on  a  guessed  or  estimated  pressure  field.  In 
the  corrector  step,  a  pressure  (or  a  pressure-correction)  equation  is  derived 
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and  solved.  Then,  the  variation  in  the  pressure  field  is  accounted  for  within  the 
momentum  equations  by  corrections  to  the  velocity  and  density  fields.  Thus, 
the  velocity,  density,  and  pressure  fields  are  driven,  iteratively,  to  better 
satisfying  the  momentum  and  continuity  equations  simultaneously  and 
convergence  is  achieved  by  repeatedly  applying  the  above-described 
procedure. 

Before  deriving  the  pressure  or  pressure  correction  equation,  the  discretized 
momentum  equations  are  first  written  in  the  following  notationally  more  suitable 
form: 


a“Up  =  +bp  -Q(VP)p.i 

NB(P) 

apVp  =  +K  -^^(VP)p.j 

NB(P) 

This  form  can  be  simplified  to 

fupl  jD[u],  0  U(VP),.11 

lv,i  1h[v]J  1  0  D[v]jl(VP),.j| 

where 


(51) 


(52) 


.  +  bp 

(VP),  =  i  JVPd£2  hW,  =  - 

ap 


(53) 


In  the  above  equations,  Q  is  the  volume  of  cell  P,  and  the  subscripts  e,  w,  n,  and 
s  refer  to  values  at  the  east,  west,  north,  and  south  faces  of  the  control  volume 
(Fig.  1 ).  Defining  the  vector  forms  of  the  above  operators  as. 


H[v]p 


H[u]p 

H[v]p 


D[u]p  0 
0  D[v]p 


(VP)p  = 


■(VP),.i' 

■(VP);’ 

.(vp),.j_ 

.(vp);_ 

(54) 


the  momentum  equations  in  vector  form  become 
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Vp-H[v],  =-D,(VP),  (55) 

Since  an  equation  for  pressure  will  be  derived  by  combining  momentum  and 
continuity,  the  discretized  form  of  the  continuity  equation  is  needed  and  is 
obtained  as  (Fig.1): 


(56) 


which  can  be  written  as 

fe^)n+A[pv.slp  =0 

^2£^)n  +  A[pU],  =0 
where 


(57) 


Uf=Vf.Sf  (58) 

For  the  calculation  of  the  mass  fluxes  across  the  control  volume  faces  (Uf)  and 
for  checking  mass  conservation,  the  values  of  the  velocity  components  are 
needed  there.  In  order  to  avoid  oscillations  which  may  result  if  a  simple  linear 
interpolation  method  is  used,  a  special  interpolation  practice  is  employed  as 
suggested  by  Rhie  [54],  Peric  [38],  and  Majumbar  [55].  The  basis  for  the 
interpolation  procedure  are  the  discretized  momentum  equations  at  the  control 
volume  centres,  as  given  by  equation  (55)  where  the  pressure  source  term  has 
been  taken  out  of  the  Q'*’  term  and  shown  explicitly.  To  evaluate  velocities  at 
the  control  volume  face  f,  terms  in  equation  (55)  are  selectively  interpolated  and 
evaluated  at  the  f  location,  according  to: 

v,-H[y],  =-D,(vp), 


(59) 
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where  the  overbar  denotes  a  linear  interpolation.  This  relation  may  be  viewed 
as  a  pseudo-momentum  equation  at  the  control  volume  face.  The  cell  face 
velocities  are  thus  made  dependent  on  the  pressure  difference  across  the  face. 
This  interpolation  practice  helps  avoiding  the  checkerboard  problems  previously 
encountered  in  collocated  variable  algorithms  [1].  For  further  details,  the  reader 
is  referred  to  Moukalled  and  Danvish  [53]. 

The  Pressure  Correction  Equation 

As  mentioned  earlier,  the  convergence  in  the  segregated  approach  is  driven  by 
the  corrector  stage  where  a  pressure  (or  a  pressure-correction)  equation  is 
solved.  Therefore,  the  first  phase  in  developing  a  segregated  solution  algorithm 
is  to  derive  such  an  equation.  The  key  step  in  the  derivation  is  to  note  that  in 
the  predictor  stage  a  guessed  or  estimated  pressure  field  from  the  previous 
iteration,  denoted  by  is  substituted  into  the  momentum  equations.  The 
resulting  velocity  field,  denoted  by  v*,  which  now  satisfies  the  momentum 
equations,  will  not  in  general  satisfy  the  continuity  equation.  Thus,  a  correction 
is  needed  in  order  to  obtain  velocity  and  pressure  fields  that  satisfy  both 
equations.  Denoting  the  pressure,  velocity,  and  density  corrections  by  P',  v'(u', 
v'),  and  p',  respectively,  the  corrected  fields  are  obtained  from: 

>  = 

.v  =  v*+v'  (u  =  u* +u',v  =  V* -i-v') 

p  =  p(">+p' 

Before  the  pressure  field  is  known,  the  velocities  obtained  from  the  solution  of 
the  momentum  equations  are  actually  u*  and  v*  rather  than  u  and  v.  Hence 
the  equations  solved  in  the  predictor  stage  are: 
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v;-H[v‘]p=-Dp(vP<"')p 

while  the  final  solution  satisfies 
Vp -H[v]p  =-D,(VP)p 

Subtracting  the  two  sets  of  equation  (62)  and  (61)  from  each  other  yields  the 

following  equation  involving  the  correction  terms: 

yi-H[v']p=-Dp(VP')p 


Moreover,  the  new  density  and  velocity  fields,  p  and  v,  will  satisfy  the  continuity 
equation  if: 

^^~^^n+A[pvs]p  =0 

Ot 

Linearizing  the  (pv)  term,  one  gets 
(p*  +  p')(v*  +  v')=  p*v*  +  p*v'  +  p'v*  +  pV' 

Substitution  of  equation  (65)  into  equation  (64)  gives 

+ a[(p*v* +p*v'  +  p'v* +pV')s]p  =0 


Rearranging,  the  following  equation  is  obtained: 

§(Pp  +Pp')+a[(p'v*  +p*v’)s]p  =^p;  -a[(p*v*)s]p  -A[(pV).S]p 
Using  equation  (63),  the  above  equation  becomes: 


|pp'+a[u-p'+p'(H[v'1-D(VP'))s]p  =-(Pf^)n-A[(pV)s]p  -a[(p'v').s] 


Finally,  substituting  density  correction  by  pressure  correction,  as  obtained  from 
the  equation  of  state,  the  new  form  of  the  pressure-correction  equation 


becomes: 
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+a[c.U-P’]^  -A|p-D(VP')sl,  =-(p'’^^P?)n-A[p-U']p  -a|p'H[v'J.s]p 

-A[(pV)5]p 

The  usual  practice  is  to  neglect  the  second  order  correction  term  p'v' .  This 
does  not  affect  neither  the  convergence  rate  (i.e.  it  is  considerably  smaller  than 
other  terms)  nor  the  final  solution,  since  at  the  state  of  convergence  the 
correction  fields  vanish.  Furthermore,  if  the  H[v']  term  in  the  above  equation  is 
retained,  there  will  result  a  pressure  correction  equation  relating  the  pressure 
correction  value  at  a  point  to  all  values  in  the  domain.  Even  though  such  an 
equation  ensures  that  the  corrected  fields  will  satisfy  the  continuity  and 
momentum  equations,  it  is  undesirable  because  it  becomes  intractable.  To 
facilitate  implementation  and  reduce  cost,  this  term  is  neglected  in  SIMPLE. 
Therefore,  the  final  form  of  the  pressure-correction  equation  is: 

^p;  +a[c.U-P'],  -A|p-D(VP')s]p  =-fc^^*£2-A[p-U‘], 

From  the  above  equation,  it  is  clear  that  the  starred  continuity  equation  appears 
as  a  source  term  in  the  pressure  correction  equation.  Moreover,  in  a  pressure- 
based  algorithm,  the  pressure-correction  equation  is  the  most  important 
equation  that  gives  the  pressure,  upon  which  all  other  variables  are  dependent. 
Therefore,  the  accuracy  of  the  predictions  depends  on  the  proper  estimation  of 
pressure  from  this  equation.  Definitely,  the  more  accurate  the  interpolated 
starred  density  (p*)  values  at  the  control  volume  faces  are,  the  more  accurate 
the  predicted  pressure  values  will  be.  The  use  of  a  central  difference  scheme 
for  the  interpolation  of  p*  leads  to  instability  at  Mach  numbers  near  or  above  1 
[12,15].  On  the  other  hand  the  use  of  a  first  order  upwind  scheme  leads  to 
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excess  diffusion  [15].  The  obvious  solution  to  the  aforementioned  problems 
would  be  to  interpolate  for  values  of  p*  at  the  control  volume  faces  in  the  same 
way  interpolation  for  other  dependent  variables  is  carried  out:  in  other  words,  to 
employ  the  bounded  HR  family  of  schemes  for  which  no  problem-dependent 
factors  are  required.  Adopting  this  strategy,  the  discretized  form  of  the  starred 
steady  continuity  equation  becomes: 


a[p-u-], =(p:ru:  4:ru>(p:ru;  4:ru; 

As  can  be  seen  from  equation  (72),  there  is  no  need  for  any  deferred  correction 
and  the  interpolated  starred  density  values  may  directly  be  used.  The  same 
procedure  is  also  adopted  for  calculating  the  density  when  computing  the  mass 
flow  rate  at  a  control  volume  face  In  the  general  conservation  equation. 

When  discretizing  the  pressure-correction  equation  (Eq.  (71)),  careful  attention 
should  be  paid  to  the  second  term  on  the  left  hand  side  that  is  similar  to  a 
convection  term  and  for  which  any  convective  scheme  may  be  used.  Since  at 
the  state  of  convergence  the  pressure-correction  field  is  zero,  the  order  of 
interpolation  scheme  is  not  important  and  does  not  affect  the  solution  accuracy. 
Therefore,  the  use  of  a  first  order  scheme  is  sufficient  for  the  purpose  of 
discretizing  this  term.  Adopting  the  UPWIND  scheme  [1],  the  discretized  form  of 
the  convection-like  term  AfCgU'P'jp  is: 

a[c.u-p'J^  =(c,Vv\  +(c,u-p1  +(c,u-p'),  +(c.u-p’), 

=  (c.)J|u;,o|p;  -||-u:.o||p;]+(c,)J|u:.o||p;  -||-u:,o||p; 

(c,)J|u;.o||p;  -|-u;,o||Pi]+(c,)j|u:,o||p;  -|-u;,o||p' 


+ 


Rearranging,  one  obtains: 
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A[c,u-pi  =  [(c.  )J|u:,o|| + (c, ),  llu;,  ,o|| + (c.  )J|u:  .o|| + (c.  )J|u;  ,0||]p;  - 

(c.)J|-u;,o||p;  -(c.)J|-u-.,o||p;  -(c,)J|-u;,o||p;  -(c,)J|-u:,o||p; 

The  term  a[p*d(VP').s]p  is  discretized  following  the  same  procedure  that  was 


used  in  discretizing  the  diffusion  flux.  Its  final  form  is  given  by: 
a[p'D(VP'):S],  =(p'D(VP')s).  +(p-D(VP')s),  +(p'D(VP')s).  +(p‘D(VP').s). 


(75) 

where  the  underlined  terms  account  for  the  non-orthogonal  factors.  They  are 
usually  neglected  since  their  contribution  is  small  in  comparison  with  other 
terms  and  vanish  when  the  grid  is  orthogonal.  However,  they  could  be 
accounted  for  by  moving  them  to  the  right  hand  side,  adding  them  to  the  source 
term,  and  modifying  the  solver  to  explicitly  update  their  values  after  a  solver 
(not  global)  iteration.  Neglecting  these  terms,  Eq.  (75)  becomes: 
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A[p-D(vp')s],  =p: '^d[v].(s;)  _p, 

o  .cl 

•^e  *^6 

•  D[ui.(s;y+D[vL(s:L  s 

Hw  C  A 

o  W 

p.Su]M±mkI(p,.p,), 

n  n 


S  .d 

s  s 


(76) 


=  rr  (PJ  -  p; ) + r;  (p;  -  p; )  +  (p;  -  p;  >  r;  (p;  -  p; ) 

= -(C  +  r:'  +  r:  +  r/'  )p;  +  rf  p'  +  r:  p;  +  r:'p;  +  r;p' 

Moreover,  the  discretized  form  of  a[p'U‘]p  is  given  by  Eq.  (72).  Substituting  the 
various  terms  in  Eq.  (71)  by  their  equivalent  expressions  as  derived  above,  the 
final  form  of  the  pressure-correction  equation  is: 


a^'p;  =a^'P'  +a^p;  +a^'p;  +a^'P'  +br 


where 

aE=rr+(c„)j|-u:,o| 

a:;=r:+(c„)j|-u:.o| 

aJ=r:+(c„)J|-u;,o|| 


a?  =C  + 


.,.,rU.,0| 


«(Cp) 


(a:)’ 

Ul 

aj-  =  a'-  +  a:;  +  a''  +  af  +  (a^  )■  +  [{C, ),  u;  +  (c,  ),  U’.  +  (c,  )„  u;  +  (c, ).  u;  ] 
b.  _k^)n-((p:)“u;  ^pTv:  4;r  u;  +{p:r  u;) 


(77) 


(78) 
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Overall  Solution  Procedure 

Knowing  the  solution  at  time  t,  the  solution  at  time  t+5t  is  obtained  as  follows: 

•  Solve  implicitly  for  u  and  v,  using  the  available  pressure  and  density  fields. 

•  Calculate  the  2)  field. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u,  v,  P  and  p. 

•  Solve  implicitly  the  energy  equation  and  update  the  density  field. 

•  Return  to  the  first  step  and  iterate  until  convergence. 

Boundary  Conditions 

The  solutions  to  the  above  system  of  equations  require  the  specification  of 
boundary  conditions  of  which  several  types  are  encountered  in  flow 
calculations,  such  as  inflow,  outflow,  and  no-flow  (impermeable  walls,  and 
symmetry  lines).  Details  regarding  the  various  types  and  their  implementation 
for  both  incompressible  and  compressible  flow  calculations  are  given  below. 

Incompressible  flows 
Inflow  boundary: 

At  an  inflow  boundary,  the  flow  enters  the  computational  domain  and  normally 
the  values  of  all  variables  are  known.  These  values  can  be  directly  substituted 
into  the  discretized  equations  for  the  boundary  control  volumes  and  thus 
nothing  special  needs  to  be  done. 


A  High  Resolution  Pressure  Based  Algorithm  for  Fluid  Flow  at  All  Speeds 


No-flow  boundary 

At  a  no-flow  boundary,  the  fluid  cannot  cross  the  boundary.  Examples  of  such 
boundaries  are  impermeable  walls  and  symmetry  lines. 

Impermeable  wall 

For  viscous  flows,  the  no-slip  condition  can  be  imposed  for  the  velocities  at  an 
impermeable  wall.  According  to  this  condition,  the  velocity  of  the  fluid  at  the  wall 
must  be  equal  to  the  velocity  of  the  wall.  For  potential  flows,  the  flow  at  a  solid 
wall  is  assumed  to  be  totally  along  the  wall.  Setting  the  normal  component  of 
velocity  to  zero  imposes  this  condition.  For  other  scalars,  either  the  scalar  value 
at  the  wall  or  its  flux  may  be  specified. 

Symmetry  line 

At  a  symmetry  line,  the  normal  component  of  velocity  and  the  normal  gradient 
of  the  parallel  component  of  velocity  are  zero. 
n.v  =  0 
(n.V|v|  =  0 

For  other  scalar  variables,  both  the  convective  and  diffusive  fluxes  are  zero. 
These  conditions  may  be  expressed  as  follows: 
n.Vcj)  =  0 

In  the  above  expressions,  n  is  the  unit  vector  normal  to  the  symmetry  line. 
Outflow  boundary 

The  boundary  at  which  the  fluid  leaves  the  computational  domain  is  called  an 
outflow  boundary.  The  most  commonly  used  practice  at  an  outflow  boundary  is 
to  assume  that  the  diffusive  flux  is  zero  and  the  total  flux  is  purely  convective  in 
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nature.  This  assumption  is  equivalent  to  setting  the  streamwise  gradients  to 
zero. 

Compressible  flows 

Inflow  boundary 

The  number  of  variables  that  can  be  specified  at  the  inflow  boundary  depends 
on  the  number  of  incoming  characteristics  at  that  surface.  For  a  two- 
dimensional  subsonic  inflow,  this  requires  the  specification  of  three  variables. 
For  a  supersonic  flow,  however,  all  variables  must  be  known  a  priori.  For  the 
case  of  subsonic  inflow  there  is  considerable  choice  as  to  which  variables  are 
specified.  For  example,  both  the  velocity  components  and  the  temperature  may 
be  specified.  In  internal  flows  (e.g.  flow  in  nozzles),  experimentation  with  this 
approach  has  shown  poor  convergence  characteristics  and  divergence  in  the 
presence  of  shock  waves  when  HR  schemes  are  used.  Better  convergence 
properties  are  achieved,  as  described  below,  by  specifying  the  total  conditions 
(i.e.  stagnation  temperature  and  pressure)  and  the  transverse  component  of 
velocity  or  the  inlet  flow  angle.  Other  variables  that  are  required  at  the  inflow 
boundary  are  obtained  by  extrapolation  from  the  interior. 

Prescribed  total  conditions  at  inlet 

The  implementation  of  this  boundary  condition,  demonstrated  at  a  boundary 
face  f,  is  accomplished  via  a  technique  that  incorporates  implicitly  the  influence 
of  pressure  on  the  velocity.  For  that  purpose,  the  total  pressure,  total 
temperature,  and  flow  angle  are  defined  as: 
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R  =P  1  + 


y-l  v.v  V-i 


T.  =T  1  + 


tan(p)=^ 

2  yRTJ  u 


The  procedure  can  easily  be  understood  by  following  the  sequence  of  events 
during  a  global  SIMPLE  iteration.  The  solution  process  starts  by  calculating  the 
velocity  components  at  the  inlet  (f)  boundary,  Uf  and  Vf,  from  equation  (82). 
Then,  the  momentum  equations  are  solved  using  Uf  and  Vf  as  specified 


boundary  velocities  for  the  current  iteration.  In  calculating  the  convective  flux  at 
the  inlet  boundary,  however,  the  velocity  used  is  not  the  one  obtained  from 
equation  (82)  but  rather,  as  described  below,  it  is  calculated  such  that  the  mass 
conservation  is  assured.  Having  solved  the  momentum  equations,  the  mass 
fluxes  are  updated  throughout  the  solution  domain,  using  the  newly  calculated 
velocity  field.  The  new  inlet  mass  fluxes  are  then  evaluated  using  the  velocity 
components  Uf  and  Vf  defined  above.  Moreover,  in  deriving  the  pressure 
correction  equation,  the  inlet  mass  flux  should  include  a  correction  term.  In 
order  to  avoid  any  possible  negative  coefficients  in  the  resulting  equation,  the 
correction  is  applied  to  the  velocity  alone,  assuming  that  the  density  is  either 
prescribed  or  treated  as  such  within  one  outer  iteration  [12,18].  Mathematically, 
the  velocity  corrections  are  expressed  via  pressure  correction  as  follows: 


Vf  =  Uf  tan(p) 


(83) 


The  coefficient 


is  obtained  from  equation  (82)  as: 


u*y(l  +  tan'(p))Pf 


y-l  u*^(l  +  tan^ 
~Y~  yRT* 


1-2y 

Y-l 


(84) 
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The  boundary  mass  flux  correction  is  expressed  as  a  function  of  the  boundary 
pressure  correction  as  follows: 


Where  P/  is  obtained  by  extrapolating  the  pressure  correction  from  interior 
nodes.  Thus,  the  coefficients  in  the  pressure  correction  equation  for  the  cells 
next  to  the  boundary  are  modified  accordingly.  Once  the  pressure  correction 
field  is  obtained  the  related  variables  (velocity  components,  pressure,  and 
density)  and  the  mass  fluxes,  which  are  used  to  calculate  the  convection  fluxes 
in  all  equations  in  the  following  iteration,  are  corrected  throughout  the  solution 
domain,  including  the  inlet  boundary.  The  newly  predicted  pressure  at  the  inlet 
boundary  is  used  in  equation  (82)  to  calculate  new  velocities  that  will  serve  as 
the  prescribed  boundary  velocities  in  the  next  iteration. 

No-flow  boundary 

The  treatment  at  a  no-flow  boundary  for  compressible  flows  is  identical  to  that 
for  incompressible  flows  discussed  above. 

Outflow  boundary 

Similar  to  the  inflow  boundary,  the  number  of  variables  that  can  be  specified  at 
an  outflow  boundary  is  equal  to  the  number  of  incoming  characteristics.  For  a 
subsonic  outflow,  this  requires  specification  of  one  boundary  condition.  The 
most  common  practice  is  to  assume  the  static  pressure  to  be  known.  All  other 
variables  are  extrapolated  from  the  interior  of  the  computational  domain. 
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Prescribed  static  pressure 

Although  the  pressure  correction  at  the  boundary  by  definition  will  be  equal  to 
zero  the  velocity  and  mass  flux  corrections  will  not.  Thus,  the  pressure 
correction  equation  needs  to  be  modified  for  the  cells  next  to  the  boundary. 

The  velocity  components  at  the  boundary  nodes  are  calculated  in  a  similar  way 
as  in  equation  (55),  namely, 
v;  =H[y],  -D,(VP), 

which  results  in  the  following  expressions  for  the  velocity  corrections: 

The  mass  flux  correction  is  calculated  as: 

(p*  V  +  p  V ),  .S,  =  p*  (-  D,  (VP')f  .S, )+  Cp^'u; 

= -  p‘  PMf(Sf)  +_D[v]f(s^)  _  p,  ^  C  ^'u; 

The  resulting  modifications  for  the  coefficients  of  interior  nodes  in  the  pressure 
correction  equation  are  easily  derived  from  the  above  expression.  The  flux 
velocity  correction U'f  at  the  boundary  is  obtained,  after  solving  the  pressure- 
correction  equation,  from  the  following  formula: 


In  the  previous  equations,  the  overbar  denotes  values  which  are  obtained  by 
extrapolation.  Moreover,  the  above  procedure  is  equally  applicable  to 
incompressible  flows  and  the  corresponding  expressions  for  the  velocities  and 
mass  fluxes  at  boundaries  can  be  derived  from  those  given  above  by  setting 
the  density  correction  to  zero. 
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For  a  supersonic  outflow,  all  variables  must  be  determined  by  extrapolation. 
The  velocity  and  mass  flux,  as  well  as  their  corrections  are  obtained  in  the 
same  way  explained  above.  However,  since  the  pressure  is  not  known,  the 
pressure  correction  at  the  boundary  node,  Pp ,  is  not  zero.  Thus,  it  should  be 
expressed  in  terms  of  the  values  at  interior  nodes  and  the  coefficients  of  the 
pressure  correction  equation  for  the  next-to-boundary  cells  must  be  modified 
accordingly. 

Results  and  discussion 

The  validity  of  the  above  described  solution  procedure  is  demonstrated  in  this 
section  by  presenting  solutions  to  the  following  four  inviscid  test  cases:  (i)  flow 
in  a  converging  diverging  nozzle;  (ii)  flow  over  a  bump;  (iii)  supersonic  flow  over 
a  step:  and  (iv)  the  unsteady  duct  filling  problem. 

Flow  in  a  converging-diverging  nozzle 

The  first  test  selected  is  a  standard  one  that  has  been  used  by  several 
researchers  for  comparison  purposes  [18,19].  The  problem  is  first  solved  using 
a  pseudo-one-dimensional  variable  area  code.  The  cross-sectional  area  of  the 
nozzle  varies  as 

where  Si=2.035  and  Sth=1  are  the  inlet  and  throat  areas,  respectively,  and  0  <  x 
<  10.  Solutions  are  obtained  over  a  wide  range  of  inlet  Mach  numbers  ranging 
from  the  incompressible  limit  (M=0.1)  to  supersonic  (M=7),  passing  through 
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transonic  and  including  strong  normal  shock  waves. 

At  all  Mach  number  values,  two  sets  of  results  are  generated  and  compared 
against  the  exact  analytical  solution.  The  first  is  obtained  using  the  third-order 
SMART  scheme  [2]  for  all  variables  except  density  (for  which  the  UPWIND  [21] 
scheme  is  used).  In  the  second  set  however,  the  SMART  scheme  is  used  for  all 
variables  including  density.  Results  displayed  in  Figs.  4,  5,  and  6  are  for  inlet 
Mach  numbers  of  0.1,  0.3,  and  7,  respectively.  In  all  three  figures,  results  are 
displayed  for  four  different  uniform  grid  networks  of  sizes  125,  79,  59,  and  39 
control  volumes.  As  expected,  results  shown  in  Fig.  4  (Min=0.1,  subsonic 
throughout)  reveal  a  decrease  in  accuracy  as  the  grid  density  decreases.  In 
addition,  for  all  grid  sizes,  the  solution  is  nearly  insensitive  to  using  a  HR 
scheme  when  interpolating  for  density.  This  is  expected,  since  for  this  inlet 
Mach  number  value,  variations  in  density  are  small  and  the  flow  can  be 
considered  to  be  nearly  incompressible.  For  Min=0.3  (Fig.  5),  the  backpressure 
is  chosen  such  that  a  supersonic  flow  is  obtained  in  the  diverging  section  (i.e. 
Mth=1,  transonic).  The  Mach  number  distributions  after  the  throat  are  depicted 
in  Fig.  5.  As  shown,  the  use  of  a  HR  scheme  for  interpolating  the  values  of 
density  at  the  control  volume  faces  improves  predictions.  This  improvement 
decreases  with  increasing  grid  density.  However,  except  for  the  value  near  the 
exit  section,  results  predicted  with  values  of  density  at  the  control  volume  faces 
calculated  using  a  HR  scheme  are  nearly  coincident  with  the  exact  solution 
even  for  the  coarsest  grid  tested  (see  Fig.  5(d)).  The  Mach  number 
distributions  depicted  in  Fig.  6  are  for  a  fully  supersonic  flow  in  the  nozzle.  The 
trend  of  results  is  similar  to  that  of  Fig.  5.  Again  important  improvements  are 
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obtained  when  using  the  SMART  scheme  for  density  interpolation. 

The  accuracy  of  the  new  technique  in  predicting  normal  shock  waves  is 
revealed  by  the  Mach  number  and  pressure  distributions  displayed  in  Figs.  7(a) 
and  7(b),  respectively.  Two  backpressure  values  that  cause  normal  shock 
waves  at  x=7  and  9  are  used.  For  each  back  pressure,  three  different  solutions 
(one  using  the  UPWIND  scheme  for  all  variables;  the  second  one  using  the 
SMART  scheme  for  all  variables;  the  third  one  using  the  SMART  scheme  for  all 
variables  except  density  for  which  the  UPWIND  scheme  is  used)  are  obtained 
and  compared  against  the  exact  solution.  All  solutions  are  obtained  by 
subdividing  the  domain  into  121  uniform  control  volumes.  As  shown,  predictions 
obtained  using  the  UPWIND  scheme  for  all  variables  are  very  smooth  but  highly 
diffusive  and  cause  a  smearing  in  the  shock  wave.  Results  obtained  using  the 
SMART  scheme  for  all  variables  except  density  are  more  accurate  than  those 
obtained  with  the  UPWIND  scheme  and  cause  less  smearing  in  the  shock 
waves.  The  best  results  are,  however,  obtained  when  employing  the  SMART 
scheme  for  all  variables  including  density.  The  plots  also  reveal  that  solutions 
obtained  using  the  SMART  scheme  show  some  oscillations  behind  the  shock. 
This  is  a  feature  of  all  HR  schemes.  The  oscillations  are  usually  centered  on  the 
accurate  solution  and  are  reduced  with  grid  refinement  in  both  wavelength  and 
amplitude. 

As  a  further  check  on  the  applicability  of  the  new  technique  in  the  subsonic, 
transonic,  and  supersonic  regimes,  results  are  generated  for  several  inlet  Mach 
number  values  0.1<Min<7  and  displayed  in  Fig.  8.  As  shown,  the  Mach  number 
and  pressure  distributions  are  in  excellent  agreement  with  the  exact  solution. 
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Two-dimensional  predictions  for  some  of  the  above-presented  cases  were 
generated  with  100x15  mesh  covering  one  half  of  the  nozzle.  The  resultant 
area-averaged  variations  of  Mach  number  are  depicted  in  Fig.  9.  Results  were 
obtained  using  the  SMART  scheme  for  all  variables  including  density.  As  for  the 
quasi-one-dimensional  predictions,  results  are  in  excellent  agreement  with  the 
exact  solutions. 

Flow  over  a  circular  arc  bump 

The  physical  situation  consists  of  a  channel  of  width  equal  to  the  length  of  the 
circular  arc  bump  and  of  total  length  equal  to  three  lengths  of  the  bump.  This 
problem  has  been  used  by  many  researchers  [18,19,23]  to  test  the  accuracy 
and  stability  of  numerical  algorithms.  Results  are  presented  for  three  different 
types  of  flow  (subsonic,  transonic,  and  supersonic).  For  subsonic  and  transonic 
calculations,  the  thickness-to-chord  ratio  is  10%  and  for  supersonic  flow 
calculations  it  is  4%.  In  all  flow  regimes,  predictions  obtained  over  a  relatively 
coarse  grid  using  the  SMART  scheme  for  all  variables  including  density  are 
compared  against  results  obtained  over  the  same  grid  using  the  SMART 
scheme  for  all  variables  except  density,  for  which  the  UPWIND  scheme  is  used. 
Due  to  the  unavailability  of  an  exact  solution  to  the  problem,  a  solution  using  a 
dense  grid  is  generated  and  treated  as  the  most  accurate  solution  against 
which  coarse  grid  results  are  compared. 

Subsonic  flow  over  a  circular  arc  bump 

With  an  inlet  Mach  number  of  0.5,  the  inviscid  flow  in  the  channel  is  fully 
subsonic  and  symmetric  across  the  middle  of  the  bump.  At  the  inlet,  the  flow  is 
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assumed  to  have  uniform  properties  and  all  variables,  except  pressure,  are 
specified.  At  the  outlet  section,  the  pressure  is  prescribed  and  all  other 
variables  are  extrapolated  from  the  interior  of  the  domain.  The  flow  tangency 
condition  is  applied  at  the  walls.  As  shown  in  Fig.  10(a),  the  physical  domain  is 
non-uniformly  decomposed  into  63x16  control  volumes.  The  dense  grid  solution 
is  obtained  over  a  mesh  of  size  252x54  control  volumes.  Isobars  displayed  in 
Fig.  10(b)  reveal  that  the  coarse  grid  solution  obtained  with  the  SMART  scheme 
for  all  variables  falls  on  top  of  the  dense  grid  solution.  The  use  of  the  upwind 
scheme  for  density  however,  lowers  the  overall  solution  accuracy.  The  same 
conclusion  can  be  drawn  when  comparing  the  Mach  number  distribution  along 
the  lower  and  upper  walls  of  the  channel.  As  seen  in  Fig.  10(c),  the  coarse  grid 
profile  obtained  using  the  SMART  scheme  for  density  is  closer  to  the  dense  grid 
profile  than  the  one  predicted  employing  the  upwind  scheme  for  density.  The 
difference  in  results  between  the  coarse  grid  solutions  is  not  large  for  this  test 
case.  This  is  expected  since  the  flow  is  subsonic  and  variations  in  density  are 
relatively  small.  Larger  differences  are  anticipated  in  the  transonic  and 
supersonic  regimes. 

Transonic  flow  over  a  circular  arc  bump 

With  the  exception  of  the  inlet  Mach  number  being  set  to  0.675,  the  grid 
distribution  and  the  implementation  of  boundary  conditions  are  identical  to 
those  described  for  subsonic  flow.  Results  are  displayed  in  Fig.  1 1  in  terms  of 
isobars  and  Mach  profiles  along  the  walls.  In  Fig.  11(a)  isobars  generated  using 
a  dense  grid  and  the  SMART  scheme  for  all  variables  are  displayed.  Fig.  11(b) 
presents  a  comparison  between  the  coarse  grid  and  dense  grid  results.  As 
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shown,  the  use  of  the  HR  SMART  scheme  for  density  greatly  improves  the 
predictions.  Isobars  generated  over  a  coarse  grid  (63x16  c.v.)  using  the 
SMART  scheme  for  all  variables  are  very  close  to  the  ones  obtained  with  a 
dense  grid  (252x54  c.v.).  This  is  in  difference  with  coarse  grid  results  obtained 
using  the  upwind  scheme  for  density  and  the  SMART  scheme  for  all  other 
variables,  which  noticeably  deviate  from  the  dense  grid  solution.  This  is  further 
apparent  in  Fig.  11(c)  where  Mach  number  profiles  along  the  lower  and  upper 
walls  are  compared.  As  shown,  the  most  accurate  coarse  grid  results  are  those 
obtained  with  the  SMART  scheme  for  all  variables  and  the  worst  ones  are 
achieved  with  the  upwind  scheme  for  all  variables.  The  maximum  Mach  number 
along  the  lower  wall  (si  .41),  predicted  with  a  dense  grid,  is  in  excellent 
agreement  with  published  values  [18,19,23].  The  use  of  a  HR  scheme  for 
density  greatly  enhances  the  solution  accuracy  with  coarse  grid  profiles 
generated  using  the  SMART  scheme  for  all  variables  being  very  close  to  the 
dense  grid  results.  By  comparing  coarse  grid  profiles  along  the  lower  wall,  the 
all-SMART  solution  is  about  11%  more  accurate  than  the  solution  obtained 
using  SMART  for  all  variables  and  upwind  for  density  and  21%  more  accurate 
than  the  highly  diffusive  all-upwind  solution. 

Supersonic  flow  over  a  circular  arc  bump 

Computations  are  presented  for  two  inlet  Mach  number  values  of  1 .4  and  1 .65. 
For  these  values  of  inlet  Mach  number  and  for  the  used  geometry,  the  flow  is 
also  supersonic  at  the  outlet.  Thus,  all  variables  at  inlet  are  prescribed,  and  at 
outlet  all  variables  are  extrapolated.  For  Min=1.4,  results  are  presented  in  Figs. 
12  and  13.  The  coarse  grid  used  is  displayed  in  Fig.  12(a).  Mach  number 
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contours  are  compared  in  Fig.  12(b).  As  before,  the  coarse  grid  all-SMART 
results  (58x18  c.v.),  being  closer  to  the  dense  grid  results  (158x78  c.v.),  are 
more  accurate  than  those  obtained  when  using  the  upwind  scheme  for  density. 
The  fine  grid  Mach  contours  are  displayed  in  Fig.  12(c).  As  depicted,  the 
reflection  and  intersection  of  the  shocks  is  very  well  resolved  without  undue 
oscillations.  The  Mach  profiles  along  the  lower  and  upper  walls,  depicted  in 
Fig.  12(d),  are  in  excellent  agreement  with  published  results  [56]  and  reveal 
good  enhancement  in  accuracy  when  using  the  SMART  scheme  for  evaluating 
interface  density  values.  The  use  of  the  upwind  scheme  to  compute  density 
deteriorates  the  solution  accuracy  even  though  a  HR  scheme  is  used  for  other 
variables.  The  all-upwind  results  are  highly  diffusive.  Finally,  results  for  this 
case  were  obtained  over  a  grid  of  90x30  nodes,  of  which  80x30  were  uniformly 
distributed  in  the  region  downstream  of  the  bump’s  leading  corner.  Resulting 
Mach  contours  are  compared  in  Fig.  13  with  four  other  solutions  [19,57,58] 
using  the  same  grid  density.  The  comparison  demonstrates  the  credibility  and 
superiority  of  the  current  solution  methodology.  The  wiggles  and  oscillations  in 
some  regions  around  the  shock  waves  in  the  published  solutions  are  not 
present  in  the  newly  predicted  one. 

For  Min=1.65,  results  are  depicted  in  Fig.  14.  The  coarse  grid  used  is  shown  in 
Fig.  14(a)  and  the  Mach  contours  are  compared  in  Fig.  14(b).  The  trend  of 
results  is  consistent  with  what  was  obtained  earlier.  Fine  grid  results  displayed 
in  Figs.  14(c)  and  14(d)  are  in  excellent  agreement  with  published  results 
[18,23].  The  Mach  contours  in  Fig.  14(c)  are  very  smooth  and  do  not  show  any 
sign  of  oscillations.  The  profiles  along  the  lower  and  upper  walls  indicate  once 
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more  that  the  use  of  a  HR  scheme  for  density  increases  the  solution  accuracy. 
Thus,  for  subsonic,  transonic,  and  supersonic  flows  the  use  of  a  HR  scheme  for 
calculating  interface  density  values  undoubtedly  increases  the  solution 
accuracy. 

Supersonic  flow  over  a  step 

The  physical  situation  and  boundary  conditions  for  the  problem  are  depicted  in 
Fig.  15(a).  The  problem  was  first  solved  using  the  upwind  scheme  and  the 
predicted  isobars  are  depicted  in  Fig.  15(b).  In  fig.  15(c),  the  isobars  reported  in 
[17]  are  presented.  As  shown,  the  current  predictions  fall  on  top  of  the  ones 
reported  by  March!  and  Maliska  [17]  eliminating  any  doubts  about  the 
correctness  of  the  implementation  of  the  solution  algorithm  and  boundary 
conditions.  The  isobars  resulting  from  a  dense  grid  solution  (23x108  c.v.)  using 
the  upwind  scheme  for  all  variables  are  presented  in  Fig.  16(a).  The 
effectiveness  of  using  a  HR  scheme  for  density  is  demonstrated  through  the 
comparison  depicted  in  Fig.  16(b).  Two  different  isobars  representing  pressure 
ratios  of  values  0.9  and  2.5  are  considered.  Solutions  obtained  over  a  course 
grid  (38x36  c.v.)  using:  (i)  the  SMART  schemes  for  all  variables,  (ii)  the  SMART 
scheme  for  all  variables  except  density  and  the  upwind  scheme  for  density,  and 
(iii)  the  upwind  scheme  for  all  variables,  are  compared  against  a  dense  grid 
solution  (238x108  c.v.)  generated  using  the  upwind  scheme  for  all  variables. 
Once  more  the  virtues  of  using  a  HR  scheme  for  density  is  obvious.  The  coarse 
grid  isobars  obtained  using  the  SMART  scheme  for  ai!  variabies,  being  nearly 
coincident  with  dense  grid  isobars,  are  remarkably  more  accurate  than  coarse 
grid  results  obtained  using  the  SMART  scheme  for  all  variables  except  density 
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and  the  upwind  scheme  for  density.  Finally,  in  order  to  show  that  the 
suggested  solution  algorithm  is  applicable  to  all  HR  schemes,  the  same 
problem  was  solved  over  a  coarse  grid  using  a  number  of  HR  schemes  and 
results  are  displayed  in  Fig.  17.  All  results  exhibit  the  same  behavior  and  are  by 
far  more  accurate  than  the  solution  obtained  with  the  upwind  scheme  reported 
in  the  literature  [17,59]. 

Ideal  unsteady  duct  filling 

Having  established  the  credibility  of  the  solution  method,  an  unsteady  process 
of  duct  filling  is  considered.  The  physical  situation  for  the  problem  consists  of  a 
duct  containing  a  gas  (y=1.4)  that  is  isentropically  expanded  from  atmospheric 
pressure.  The  duct  is  considered  to  be  frictionless,  adiabatic,  and  of  constant 
cross-section.  Moreover,  it  is  assumed  that  the  duct  is  opened  instantaneously 
to  the  surrounding  atmosphere,  inflow  is  isentropic,  and  in  the  fully  open  state 
the  effective  flow  area  at  the  duct  end  is  equal  to  the  duct  cross-sectional  area. 
The  unsteady  one-dimensional  duct  filling  process  is  solved  using  a  two- 
dimensional  code  over  a  uniform  grid  of  density  299x3  control  volumes,  a  time 
step  of  value  10"'^,  and  the  SMART  scheme  for  all  variables. 

The  problem  is  solved  for  a  surrounding  to  duct  pressure  ratio  of  2.45  and 
generated  results  are  displayed  in  Fig.  18.  Due  to  the  lower  pressure  of  the  gas 
contained  in  the  duct,  when  the  duct  end  is  suddenly  opened,  a  compression 
wave  is  established  instantly  as  a  shock  wave.  The  wave  diagram  for  the 
process  is  shown  in  Fig.  18(a).  The  shock  wave  moves  in  the  duct  until  the 
closed  end  is  reached.  On  reaching  the  closed  end,  the  compression  wave  is 
reflected  and  the  duct  filling  process  continues  until  the  reflected  shock  wave  is 
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at  the  open  end.  Beyond  that,  duct  emtying  starts  and  computations  were 
stopped  at  that  moment  in  time.  In  addition,  the  path  of  the  first  particle  to  enter 
the  duct  is  shown  in  the  figure.  This  was  computed  by  storing  the  duct  velocities 
at  ali  time  steps  and  then  integrating  in  time  to  locate  the  position  of  the  particle. 
Results  depicted  in  Fig.  18(a)  were  compared  against  similar  ones  reported  by 
Azoury  [60]  using  a  graphical  method.  The  two  sets  of  results  were  found  to  be 
in  excellent  agreement  with  the  ones  computed  here  falling  right  on  top  of  those 
reported. 

The  variation  of  Mach  number  with  time  at  the  open  end  of  the  duct  is  diplayed 
in  Fig.  18(b).  With  the  exception  of  the  slight  overshoot  at  the  beginning  of  the 
computations,  the  Mach  number  remains  constant  throughout  the  filling  process 
and  it  instantaneously  decreases  to  zero  at  the  time  when  the  reflected  shock 
wave  reaches  the  open  end  of  the  duct.  When  using  the  same  reference 
quantities,  the  analytical  solution  to  the  problem  reported  in  [60]  predicts  a 
constant  Mach  number  of  value  0.4391  which  is  0.21%  different  than  the  one 
obtained  here.  Moreover,  the  instantaneous  decrease  of  Mach  number  to  zero 
is  well  predicted  by  the  method.  Finally,  the  increase  in  mass  within  the  duct  is 
presented  in  Fig.  18(c).  As  expected,  due  to  the  constant  value  of  the  inlet 
Mach  number  the  mass  increases  linearly  with  time. 

Concluding  Remarks 

A  new  collocated  high-resolution  pressure-based  algorithm  for  the  solution  of 
fluid  flow  at  all  speeds  was  formulated.  The  new  features  in  the  algorithm  are 
the  use  of  a  HR  scheme  in  calculating  the  density  values  at  the  control  volume 
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faces  and  the  use  of  the  NVSF  methodology  for  bounding  the  convection 
fluxes.  The  method  was  tested  by  solving  four  problems  representing  flow  in  a 
converging-diverging  nozzle,  flow  over  a  bump,  flow  over  an  obstacle,  and 
unsteady  duct  filling.  Mach  number  values  spanning  the  entire  subsonic  to 
supersonic  spectrum,  including  transonic  flows  with  strong  normal  shock 
waves,  were  considered.  In  all  cases,  results  obtained  were  very  promising  and 
revealed  good  enhancement  in  accuracy  at  high  Mach  number  values  when 
calculating  interface  density  values  using  a  High-Resolution  scheme. 
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Figure  Captions 

Fig.  1  Control  Volume. 

Fig.  2  Typical  control  volume  faces  and  geometric  nomenclature. 

Fig.  3  (a)  Control  volume  nodes;  (b)  Normalization:  (c)  CBC  Criterion. 

Fig.  4  Comparison  of  Mach  number  variation  for  an  inlet  Mach  number  of 
0.1. 

Fig.  5  Comparison  of  Mach  number  variation  for  an  inlet  Mach  number  of 
0.3. 

Fig.  6  Comparison  of  Mach  number  variation  for  an  inlet  Mach  number  of  7. 

Fig.  7  Transonic  inviscid  flow  in  a  converging-diverging  nozzle:  (a)  Mach 
number  distributions,  (b)  pressure  distributions. 

Fig.  8  Comparison  of  (a)  Mach  number  and  (b)  pressure  distributions  for 
one-  dimensional  inviscid  nozzle  flow. 

Fig.  9  Comparison  of  distributions  of  (a)  area-averaged  Mach  number  and 
(b)  pressure  for  inviscid  nozzle  flow  from  2-D  solution. 

Fig.  10  Subsonic  flow  over  a  10%  circular  bump;  (a)  coarse  grid  used,  (b) 
isobars,  and  (c)  profiles  along  the  walls. 

Fig.  11  Transonic  flow  over  a  10%  circular  bump;  (a)  Isobars  using  a  dense 
grid,  (b)  isobars  using  various  schemes,  and  (c)  profiles  along  the 


A  High  Resolution  Pressure  Based  Algorithm  for  Fluid  Flow  at  All  Speeds 


53 


walls. 

Fig.  12  Supersonic  flow  over  a  4%  circular  bump  (Min=1.4);  (a)  coarse  grid 
used,  (b)  Mach  number  contours  using  various  schemes,  (c)  Mach 
number  contours  using  a  dense  grid,  and  (d)  profiles  along  the  walls. 

Fig.  13  Supersonic  inviscid  flow  over  4%  bump  (Min=1.4):  Mach-number 
contours. 

Fig.  14  Supersonic  flow  over  a  4%  circular  bump  (Min=1.65);  (a)  coarse  grid 
used,  (b)  Mach  number  contours  using  various  schemes,  (c)  Mach 
number  contours  using  a  dense  grid,  and  (d)  profiles  along  the  walls. 

Fig.  15  Supersonic  flow  over  an  obstacle;  (a)  Physical  situation,  (b)  Isobars 
using  the  upwind  scheme  (40X38  grid  points),  and  (c)  results  obtained 
by  Marchi  and  Maliska  using  (44x36  grid  points). 

Fig.  16  Supersonic  flow  over  an  obstacle:  (a)  Isobars  generated  using  a 
dense  grid;  (b)  Isobars  generated  using  different  schemes  for 
supersonic  flow  over  an  obstacle. 

Fig.  17  Isobars  for  the  flow  over  an  obstacle  obtained  from  different  HR 
schemes. 

Fig.  18  (a)  Wave  diagram  for  optimum  duct-filling  process;  (b)  Mach  number 

distribution  at  inflow;  (c)  Variation  of  mass  with  time. 
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Fig.  2  Typical  control  volume  faces  and  geometric  nomenclature. 


Fig.  4  Comparison  of  Mach  number  variation  for  an  inlet  Mach  number  of  0. 1 . 
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Fig.  7  Transonic  inviscid  flow  in  Laval  nozzle:  (a)  Mach  number  distributions,  (b)  pressure 


distributions. 


Fig.  10  Subsonic  flow  over  a  10%  circular  bump;  (a)  coarse  grid  used,  (b)  isobars,  and  (c)  profiles 


along  the  walls. 


Fig.  12  Supersonic  flow  over  a  4%  circular  bump  {Min=1.4);  (a)  coarse  grid  used,  (b)  Mach  number 
contours  using  various  schemes,  (c)  Mach  number  contours  using  a  dense  grid,  and  (d) 
profiles  along  the  walls. 


Fig.  13  Supersonic  inviscid  flow  over  4%  bump  (Min=1 .4):  Mach-number  contours. 
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(d) 

Fig.  14  Supersonic  flow  over  a  4%  circular  bump  (Min=1.65);  (a)  coarse  grid  used,  (b)  Mach 
number  contours  using  various  schemes,  (c)  Mach  number  contours  using  a  dense  grid, 
and  (d)  profiles  along  the  walls. 


(C) 

Fig.  15  Supersonic  flow  over  an  obstacle:  (a)  Physical  situation,  (b)  Isobars  using  the  upwind 
scheme  (40X38  grid  points),  and  (c)  results  obtained  by  Marchi  and  Maliska  using  (44x36 
grid  points). 


Fig.  17  Isobars  for  the  flow  over  an  obstacle  obtained  from  different  HR  schemes. 


Fig.  18  (a)  Wave  diagram  for  optimum  duct-filling  process,  (b)  Mach  number  distribution  at  inflow;  (c)  Variation  of  mass  with  time. 


